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This is a status report on our endeavor to reveal the mechanism of core-collapse super- 
novae (CCSNe) by large-scale numerical simulations. Multi-dimensionality of the supernova 
engine, general relativistic magnetohydrodynamics, energy and lepton number transport by 
neutrinos emitted from the forming neutron star as well as nuclear interactions there, are all 
believed to play crucial roles in repelling infalling matter and producing energetic explosions. 
These ingredients are nonlinearly coupled with one another in the dynamics of core-collapse, 
bounce, and shock expansion. Serious quantitative studies of CCSNe hence make extensive 
numerical computations mandatory. Since neutrinos are neither in thermal nor in chemical 
equilibrium in general, their distributions in the phase space should be computed. This is a 
six dimensional (6D) neutrino transport problem and quite a challenge even for those with 
an access to the most advanced numerical resources such as the "K computer" . To tackle this 
problem, we have embarked on multi-front efforts. In particular we report in this paper our 
recent progresses in the treatments of multi-dimensional (multi-D) radiation- hydrodynamics. 
We are currently proceeding on two different paths to the ultimate goal; in one approach 
we employ an approximate but highly efficient scheme for neutrino transport and treat 3D 
hydrodynamics and/or general relativity rigorously; some neutrino-driven explosions will be 
presented and comparisons will be made between 2D and 3D models quantitatively; in the 
second approach, on the other hand, exact but so far Newtonian Boltzmann equations are 
solved in two and three spatial dimensions; we will show some demonstrative test simula- 
tions. We will also address the perspectives of exa-scale computations on the next generation 
supercomputers. 

§1. Explorations of CCSNe by multi-D simulations 

From the very beginning,^) CCSN research has been one of the greatest chal- 
lenges in computational astrophysics. This is simply due to the following facts (see, 
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e.g., 2), 3) and references therein): 3D macroscopic magnetohydrodynamics on stel- 
lar scales is largely dictated by weak/strong interactions on femto-meter scales; 
not only dynamical time scales of milliseconds but neutrino-diffusion time scales 
of > 100ms become also important in the dynamics of massive star cores; compact- 
ness of proto-neutron stars (PNSs) and high velocities of collapsing material make 
it simply impossible to neglect special/general relativity; (magneto)hydrodynamical 
instabilities render the core dynamics intrinsically non-spherical; phase transitions 
in dense hadronic matter and/or neutrino oscillations may have a critical impact. 
These diverse physical elements involved in CCSNe, although they are the reasons 
why many researchers from different disciplines have been attracted so much, are 
the very sources of difficulties in numerical simulations. 

Energetics also suggests the subtlety in the explosion dynamics. Recently ad- 
vancing ability of the HST has enabled the direct observation of the progenitors of 
nearby CCSNe from pre-explosion images (see^^ for review). The accumulating ob- 
servations suggest that the majority of CCSNe (the so-called type Il-plateau (II-P)) 
comes predominantly from stars in the range about of 8 - 16 Mq.^^ A generic ex- 
plosion energy for the progenitors in the mass range is roughly on the order of 10^^ 
erg,^) which is much smaller than the available energy, ~ 10^^ erg, i.e., the gravita- 
tional binding energy of PNSs. This large amount of energy is stored as the internal 
energy of PNSs and tapped slowly by neutrinos. After the successful detections 
of neutrinos from SN1987A followed by detailed analyses of these events (e.g., 6), 
and see also 7) for a recent review), CCSN researchers have a confidence that the 
explosion energy is indeed supplied by this reservoir. Present state-of-the-art super- 
nova simulations require more than 10^*^ operations and typically more than several 
months for a single model on currently available best supercomputing platforms. The 
energy conservation should be satisfied within an error of at least less than 1% 
(= lO^^erg/lO^^ erg) in such very long-term simulations to obtain reliable results. It 
may not be difficult then to imagine how demanding the numerical simulations of 
CCSNe are. 

It is interesting that after -|-45 years of intensive and extensive theoretical studies 
we are still working on the same scenario that Colgate and White (1966) envisaged 
in their seminal paper, in which they reported the first CCSN simulation. The ab- 
stract of their paper ended with the sentence, "The energy release corresponds to the 
change in gravitational potential of the unstable imploding core; the transfer of en- 
ergy takes place by the emission and deposition of neutrinos". As is now well known, 
this is exactly the essence of the so-called neutrino-heating mechanism (e.g., 8) for 
review). Although they proposed originally that the mechanism works promptly af- 
ter bounce, Bethe and Wilson later amended it^)'^^) to the currently prevailing form, 
in which a bounce-generated shock wave is stalled first, but revived later by the de- 
position of neutrinos and an explosion follows in several hundred milliseconds after 
bounce. The scenario was investigated intensively for the first 5 years in the new 
millennium under the assumption of spherically symmetry but with the full Boltz- 
mann treatment of neutrino transfer. General relativity and/or various 
neutrino reactions, some of which are supposed to be of minor importance, 
were also implemented. 
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These explorations made clear, however, that the neutrino-heating mechanism 
fails to produce explosions in ID spherical symmetry except for super- AGB stars at 
the low-mass end.^^-* Not deterred by this failure, researchers changed gear to multi- 
D modeling. By this time there had already been mounting observational evidence 
that supernova explosions are indeed aspherical in general (see, e.g., 19)-21) and 
references therein). Numerical experiments also suggested that breaking spherical 
symmetry holds a key to success of the neutrino-heating mechanism; convective 
motions (e.g., 22)-26)) and/or the so-called standing accretion shock instability, 
or SASI, (e.g., 27)-35) and references therein) help the onset of neutrino-driven 
explosions. 

In the following years, we have indeed witnessed some exploding models by the 
neutrino-heating mechanism in axisymmetric 2D simulations (see, e.g., table 1 in 3)). 
Employing one of the best approximations for 2D neutrino transfer, Buras et al.^^^ 
reported explosions firstly for a non-rotating low-mass (II.2M0) progenitor^^) and 
then for a 15AfQ progcnitor"^^) with a moderate rotation being imposed. Imple- 
menting a multi-group flux-limited diffusion algorithm to their CHIMERA code in a 
ray-by-ray manner, on the other hand, Bruenn et al.^°) obtained explosions for non- 
rotating progenitors^^) in the mass range from 12Mq to 25Mq. Implementing the 
ray-by-ray isotropic diffusion source approximation (IDSA)^^-* in the ZEUS code with 
a reduced set of weak interactions, Suwa et al.^^^ pointed out that a rapidly rotating 

13M0 progenitor produced a stronger explosion than the non-rotating counterpart 
did.43) 

Accompanying these successes are new questions, however. In addition to the 
apparent contradictions among the groups (see Table 1 in 3)), the models mentioned 
above produced generically under-energetic explosions at the end of simulations, with 
the diagnostic explosion energy being smaller by one or two orders of magnitude than 
the canonical kinetic energy of supernova ejecta {r^ lO^^erg). Hence, it is a legitimate 
concern whether we can obtain energetic explosions comparable to observations by 
the neutrino-heating mechanism with appropriate nucleosynthetic yield, which is 
one of the most important obscrvablcs.'^^) In the above-mentioned computations, 
the softest version of Lattimer &: Swesty's (LS) equation of state (EOS) '^^^ with an 
incompressibility at nuclear density, K, of 180 MeV, was commonly employed. In 
addition to the fact that recent experiments suggest a stiffer EOS with K = 240 ±20 
MeV,^^) it is now thought to be a serious flaw that the LS180 EOS cannot support 
a 2Mq cold neutron star that is certainly existent in the universe^ *^ . Employing a 
stiffer EOS with K = 263MeV based on the Hartree-Fock approximation,^'^) Marek 
et al.^^) found no explosion for the same progenitor model, whereas they indeed 
obtained an explosion for the Shen's EOS that is even stiffer with K = 281MeV.^^) 
Suwa et al. also pointed out that not only the incompressibility but the symmetry 
energy also matters for the success of neutrino-driven explosions. Impacts of more 
detailed properties of nuclear EOSs (such as the density dependence of symmetry 
energy and the skewness of compressibility^^-*' ^^)) on the multi-D neutrino- heating 
mechanism are remaining to be understood. 



*) The maximum mass for the LS180 EOS is about I.8M0 (see, e.g., 48), 49)). 
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The paramount interest of supernova researchers at present is 3D effects, how- 
ever. We know in fact that SASI is qualitatively different between 2D and 3D ^^)'^^) 
and it is naturally expected that this may have some consequences to success of the 
neutrino-heating mechanism in 3D. So far experimental simulations are contradict- 
ing each other: Nordhaus et al.^^^ claimed that 3D dynamics will make shock revival 
easier than 2D. The assertion was challenged later by Hanke et al.,^^) who found 
little difference in the critical neutrino luminosity for shock revival between their 2D 
and 3D simulations. In both of the experimental computations, the neutrino transfer 
was not solved and the controversy will be settled only by detailed self-consistent 3D 
simulations. In the next section, we present our recent findings that illuminate 3D 
effects on the neutrino-driven mechanism. The former part of §2 is devoted to our 
3D Newtonian hydrodynamical simulations with spectral neutrino transport whereas 
in the latter half we show the latest results of our fully general relativistic (GR) 3D 
simulations that employ a more approximate neutrino transport scheme. 

Although we will focus on the neutrino-heating mechanism in this paper, it 
should be mentioned here that there are some other viable mechanisms. In the so- 
called acoustic mechanism,^^) oscillations of PNSs are supposed to produce pressure 
perturbations and send acoustic powers to the stalled shock until it revives and 
produces an explosion. The mechanism will then be particularly important in the 
later postbounce phase, when the neutrino luminosity has already declined and the 
neutrino-heating mechanism has no chance of success. The merit of this mechanism 
is that matter accretion, the source of acoustic powers, will last long, possibly until 
the shock is revived. The scenario was challenged by Quataert et al.,^^^ however, 
who demonstrated that the amplitudes of g-mode oscillations of PNSs will not be so 
large (see also 60)). Although the additional energy input by acoustic waves is very 
appealing, it still remains an issue under vivid debates and has yet to be confirmed 
by other groups. '^^^ The magnetohydrodynamical (MHD) mechanism taps rotational 
energies of stellar cores (e.g., 61)-73), 75). See also 2) for collective references). 
Magnetic fields are expected to be amplified spontaneously by the magneto-rotational 
instability (MRI) even if they are tiny prior to collapse. This mechanism requires 
rapid rotation of stellar cores at the onset of core collapse. Recent stellar evolution 
models predict that such a condition can be realized only in the special case that 
experiences the so-called chemically homogeneous evolution^^)'^^) and applies to just 
a fraction (~ 1%) of massive stars. Further investigations arc currently hampered 
by the fact that high numerical resolutions are required to accurately compute the 
growth of the MRI.^^)'^°) We finally list other possibilities proposed in the literature: 
exotic physics in the proto-ncutron star,^^)~^^) viscous heating by MRI,^^)'^^) and 
dissipation of Alfven waves. '^^^ 

Improvements of input physics are another important ingredient in the CCSN 
simulations. One of the authors provided an EOS table to the society,^^)"^^) which 
is referred to as Shen's EOS and is based on the relativistic mean field theory and 
Thomas-Fermi approximation and is now a standard choice for core-collapse sim- 
ulations. Besides this and another representative EOS by Lattimer & Swesty,"^^-* 
new sets of EOS's have been reported recently.^°)"^^) We have joined this effort to 
expand the inventory of EOS's both above and below nuclear density for the last 
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few years. Above saturation density, we have later included hyperons in the same 
framework.^^) We have also combined the original table with an EOS for quark 
matter, adopting the MIT bag model and the Gibbs condition for the first order 
phase transition. These modifications manifest themselves at high densities and 
will be more important for CCSNe with black hole formations than those with neu- 
tron star formations. It is emphasized that the EOS at high densities can be probed 
by neutrino and/or gravitational wave signals from core-collapse events. ^^^'^^^"^'''^ 
The EOS below nuclear saturation density is no less important. In the conven- 
tional EOSs,^5),87)-89) ^i^g so-called single nucleus approximation was employed, in 
which thermally populated heavy nuclei are represented by a single, supposedly the 
most abundant nucleus. Combining experimental nuclear mass data and a mass 
formula, we have solved Saha-likc equations to obtain populations of various nuclei 
in constructing the EOS.*^^^ In doing so, not only the excluded- volume effect but 
the emergence of pasta phases as well as the modifications of bulk. Coulomb and 
surface energies by surrounding nuclcons and nuclei are also taken into account phe- 
nomenologically. We are currently preparing an EOS table including the electron 
capture rates according to the obtained populations.^^) Detailed comparisons with 
other tables^°)'^^) will be also published soon. 

The ultimate goal of CCSN simulations is 3D neutrino-radiation- (magneto-) hy- 
drodynamics in full GR, in which the exact Boltzmann equations are solved and all 
the relevant weak interactions are included with sufficient realism. It is worth men- 
tioning that more massive stars than the SN II-p mentioned earlier are expected to 
lose much of their mass and explode as hydrogen-stripped SNc (Ib/c and lib). Among 
them, the type Ic-BL SNe, which are associated with long gamma-ray bursts,^°°) aU 
show much broader lines than SNe Ic. Due to the large kinetic energies of 2 — 5 x 10^^ 
erg, they have been referred to as "hypernovae" (e.g., Ref.^'^^) for a recent review), 
the central engine of which is likely to be associated with massive stellar core-collapse 
with black-hole formation (e.g., coUapsar, see collective references in Refs.^^^^'^^^^). 
In such a system, the implementation of full GR in numerical simulations is essen- 
tial. In one of our approaches mentioned above, 3D hydrodynamics in full GR is 
first addressed in §2 with neutrino transport being approximated one way or an- 
other. On the other front, we are pursuing the Boltzmann transport in 3D first. We 
present the current status of our efforts along this path in §3. Our method is based 
on an implicit finite diff'erencing of the Boltzmann equations and the inversion of 
large matrices in a very efficient way is one of the major challenges. As discussed in 
the final section, the 3D version of the code may work only on the next generation, 
exa-scale platforms. In this sense, what we provide in the following sections is just 
a snapshot of a long, on-going documentary film that will record our struggles to 
make the "dream simulation" come true. 
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§2. 3D hydrodynamical simulations with approximate neutrino 

transport 

2.1. 3D Newtonian simulations with spectral neutrino transport 

It is generally very computationally expensive to solve neutrino transport in 3D 
and a light-bulb scheme^*^^) has been widely used so far, in which neutrino heating 
and cooling are treated in a parametric manner to trigger 3D explosions. Using 
this prescription, Nordhaus et al.^^^ was the first to argue that the critical neutrino 
luminosity for producing neutrino-driven explosions becomes significantly smaller in 
3D than in 2D (see, however, 57)). They employed the CASTRO code with an 
adaptive mesh refinement technique, by which unprecedentedly high resolution 3D 
calculations were made possible. 

Since the light-bulb scheme can capture fundamental properties of neutrino- 
driven explosions (albeit on a qualitative basis), it is one of the most prevailing 
approximations adopted in recent 3D models (e.g., 32), 33), 103)). A number of im- 
portant findings have been reported recently in these simulations, such as a potential 
role of non-axisymmetric SASI flows in generating spins (see 103), 104) as well as 
55), 105)) and magnetic fields^^^^ of pulsars, stochastic nature of gravitational-wave 
(e.g., 3), 107), 108)) and neutrino emission (sec 109) for recent review). 

To go up the ladders beyond the light-bulb scheme, we studied 3D effects on 
the supernova mechanism by performing the first 3D, multi-energy-group, radiation- 
hydrodynamical core-collapse simulations. -"^^^^ For the spectral transport, the IDSA 
scheme is implemented, which can be done rather in a straightforward manner by 
extending our 2D modules42),52) rj.^^^ 

can be made possible because we ap- 
ply the so-called ray-by-ray approach (e.g., 36)) in which the neutrino transport is 
solved along a given radial direction by assuming the angle-averaged matter distribu- 
tion in a spherically symmetric manner. From a technical point of view, it is worth 
mentioning that the ray-by-ray treatment is highly efficient in paralellization*^ on 
present supercomputers, most of which employ the message-passing-interface (MPI) 
routines. The IDSA scheme splits the neutrino distribution into two components, 
each of which is solved with different numerical techniques (sec 41) for more de- 
tails). A drawback in the current version of the IDSA scheme is that heavy lepton 
neutrinos {vx, i-e., i^^, Vt and their anti-particles) as well as the energy-coupling 
weak interactions have yet to be implemented. The approximation level of the IDSA 
scheme is basically the same as the one of the Multi-Group Flux-Limited Diffusion 
MGFLD scheme. The main advantage of the IDSA scheme is that the fluxes in the 
transparent region can be determined by the non-local distribution of sources rather 
than the gradient of the local intensity like in MGFLD. In the following, we briefly 
summarize the main results on our 3D simulations, in which we obtained a first 3D 
explosion for an 11.2 Mq star.^^^ 



along each radial ray 
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2.1.1. 3D dynamics from core-collapse through postbounce turbulence till explosion 

Figure 1 shows three snapshots, which are helpful to characterize hydrodynamic 
features in 3D simulations. The top panel corresponds to t = 15 ms after bounce, 
showing that the bounce shock stalls (indicated by inward arrows in the top right 
panel) at a radius of 150 km. Note that the colors of velocity arrows are chosen so that 
they would change from yellow to red as the absolute values become larger. Looking 
carefully at the top right panel, we find that matter flows supersonically (indicated 
by reddish arrows) into the standing shock (the central transparent sphere), and 
then advects subsonically (indicated by yellowish arrows) onto the proto-neutron 
star (PNS, the central bluish region in the top left panel). For the non-rotating 
progenitor, the dynamics of collapsing iron core proceeds perfectly spherically till 
the stall of the bounce shock. This is the reason why multi-D effects are invisible 
in the entropy (top left panel) and density (top right panel) distributions right after 
bounce. 

The middle panels show the epoch (t = 65 ms) when the neutrino-driven convec- 
tion is already active. From the right panel, turbulent motions can be seen (arrows in 
random directions) inside the standing shock, which is indicated by the boundary be- 
tween red and yellow arrows. The entropy behind the standing shock becomes high 
by the neutrino- heating (reddish regions in the left panel). The size of neutrino- 
heated hot bubbles becomes larger in a non-axisymmetric way later on, which is 
indicated by smaller structures encompassed by the stalled shock (i.e., inside the 
central greenish sphere in the left panel). 

The bottom panels {t = 125 ms) show the epoch when the revived shock is 
expanding aspherically, which is indicated by the outgoing yellowish arrows in the 
right panel. The asphericity of the expanding shocks could be more clearly visible by 
the sidewall panels. Prom the entropy distribution (left panel), the expanding shock 
is shown to touch a radius of 500 km (the projected back bottom panel). Inside the 
expanding shock (enclosed by the greenish membrane in the left panel), the bumpy 
structures of the hot bubbles are seen. In contrast to these smaller asphericities, the 
deformation of the shock surface is mild. This is a consequence of SASI, leading to 
the shock deformation dominated by low spherical-harmonics modes {i = 1,2). 

Figure 2 shows the net neutrino heating rate (left panel) and Trcs/Thcat: the 
ratio of the residency time scale to the neutrino-heating time scale (right panel) 
for the snapshot oi t = 125 ms in Fig. 1. This time scale ratio is known to be 
a useful quantity to diagnose the success (Ti-es/Tiicat ^ 1) i-C-, the neutrino-heating 
timcscalc is shorter than the dwell time scale of material in the gain region*^ or failure 
(Tres/^hcat ^ 1) of the neutrino-driven explosion (e.g., 34), 84), 111), 112)). The left 
panel of Fig. 2 shows that there forms the so-called gain region, in which the neutrino 
heating dominates over cooling (seen as reddish regions in the wall panels). As seen 
in the right panel, the time scale ratio reaches ~ 2 in the gain region, the evidence 
that the shock-revival is driven by the neutrino-heating mechanism. 



where the neutrino heating dominates over the neutrino coohng (see Ref^^^^ for more detail). 
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Fig. 1. 3D plots of entropy per baryon (left panels) and logarithmic density (right panels, in unit 
of g/cm^) for three snapshots (top; f = 15 ms, middle; f = 65 ms, and bottom; t — 125 ms 
after bounce {t = 0)) during the evolution of a (non-rotating) exploding 3D model of a 11.2 
Mq star (figures taken from 110)). In the right panels, velocities are indicated by arrows. The 
color contours in the a; = (back right), j/ = (back bottom), and z = (back left) planes are 
projected on the sidewalls of the graphs. For each snapshot, the linear scale is shown along the 
axis in unit of km (figures taken from 110), reproduced by permissions of the AAS). 



Supernova as Supercomputing Science 



9 




Fig. 2. Same as Figure 1 but for the net neutrino heating rate (left panel, logarithmic in unit 
of erg/cm'^/s) and Trcs/Thcat (right panel), which is the ratio of the residency to the neutrino 
heating time scales (see the text for details). The gain region (colored in red in left panel), in 
which neutrino heating dominates over cooling, is shown to be formed. The right panel shows 
that the condition of Trcs/Thcat 1 is satisfied behind the globally aspherical shock, which is a 
characteristics of SASI (figures taken from 110), reproduced by permissions of the AAS). 

2.1.2. 2D vs. 3D: which is more advantageous for the neutrino-driven explosion? 

The left panel of Fig. 3 shows the comparison of mass-shell trajectories between 
the 3D (red lines) and corresponding ID models (green line). At around 300ms after 
bounce, the average shock radius for the 3D model exceeds 1000km. On the other 
hand, no explosion is obtained for the ID model. The right panel of Fig. 3 shows a 
comparison of the average shock radii as a function of the postbounce time. In the 
2D model, the shock expands rather continuously after bounce. These trends in the 
ID and 2D models are qualitatively consistent with those found in 36)*^ 

Comparing the shock evolutions between the 2D (green line in the right panel of 
Fig. 3) and 3D (red line) models, we find that the shock expands much faster in 2D. 
The pink line labeled by "3D low" is the result of the low resolution 3D simulation, 
in which the azimuthal grid number is reduced to half the number for the standard 
model. Note that the 3D computational grid consists of 300 logarithmically spaced 
radial zones and 64 polar (6) and 32 azimuthal {(p) uniform mesh points to cover the 
entire sphere with a radius of 5000km. Compared with the standard 3D model (red 
line), the shock expansion is less energetic for the low resolution model (later than 
~ 150ms). These results indicate that a successful explosion is most easily obtained 
in 2D and hampered by low resolutions. At first glance, this may be at odds with 

*•* The reason why the shock of our 2D model expands on average much faster than theirs might 
be our neglect of general relativity, inelastic neutrino-electron scattering and cooling by heavy- 
lepton neutrinos. All of these important missing ingredients in our 3D simulations could give a 
more optimistic condition for explosions. Apparently these ingredients should be appropriately 
implemented, which we hope to be tractable in the next-generation 3D simulations. 
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Fig. 3. Left panel: the time evolution of 3D model, visualized by mass shell trajectories in thin 
gray lines. Thick red lines show the positions of shock: the maximum (top), average (middle) 
and minimum (bottom) radii. The green line represents the shock position for the ID model. 
The numbers "1.30" and "1.40" indicate the masses in unit of Mq that are enclosed inside each 
mass-shell. Right panel: the evolution of average shock radii for the 2D (green line) and 3D 
(red line) models. The pink line presents the low resolution 3D model, in which the azimuthal 
grid number is reduced to half the number for the standard model (see the body). Figures are 
taken from 110) (reproduced by permissions of the AAS). 



the results obtained in the parametric 3D explosion models (e.g., 56)), in which the 
authors claimed that explosions would be easier in 3D than in 2D. The reason for 
the discrepancy is summarized shortly. 

Fig. 4 compares the blast morphologies in the 3D (left panel) and 2D (right) 
models. In the former, non-axisymmetric structures are clearly seen. Performing a 
tracer-particle analysis, we find that the maximum residency time is longer in 3D 
than in 2D owing to the non-axisymmetric motions (see Fig. 5). This is one of ad- 
vantageous aspects of 3D models to obtain the neutrino-driven explosions. Another 
merit in 3D is that convective matter motions below the gain radius is much more 
violent than in 2D, which enhances the neutrino luminosity in 3D (see 110) for more 
details). The negative point, on the other hand, is lower energies of emitted neu- 
trinos owing to the enhanced cooling. The competition of these effects eventually 
leads to a longer neutrino-heating time scale in our 3D models with an outcome 
of a smaller net-heating rate compared with the corresponding 2D model (Fig. 6). 
Note here that the IDSA scheme, with which the feedback from the mass accretion 
to the neutrino luminosity is automatically and self-consistently incorporated unlike 
the light-bulb approximation that assumes a constant luminosity, is quite efficient 
and a good choice for the first-generation 3D simulations. 

Although it is encouraging that the shock expansion becomes more energetic with 
better resolution (recall that the explosions obtained so far are all under-energetic^^ 
and the present model is no exception), this implies that a systematic and time con- 
suming convergence test is required to draw a robust conclusion (e.g. 57)). More 
advanced treatments of neutrino transport as well as of gravity will be also needed, 
which will probably be a subject done on forthcoming petaflops-class supercomput- 
ers. 
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Fig. 4. Blast morphologies at t = 178ms after bounce in our 3D (left) and 2D (right) models 
presented as the volume rendering of entropy. The progenitor is a 11.2M0 star.^*'') The polar 
axis is tilted by about 7r/4 in both panels (figures taken from 110), reproduced by permissions 
of the AAS). 




Fig. 5. Left panel: streamlines of selected tracer particles advecting through the shock to the PNS 
in the 3D model. The observer is located on the polar axis. Several surfaces of constant entropy 
are also displayed. An outer greenish one marks the shock surface whereas the PNS is presented 
as a central sphere. The linear scale is given at the bottom right corner. Right panel: the 
number of tracer-particles traversing the gain region as a function of their individual residency 
time for the 2D and 3D models given at t = 100ms after bounce. The maximum residency time 
in the 3D model is longer than that in the corresponding 2D model (figures taken from 110), 
reproduced by permissions of the AAS). 
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Fig. 6. Time evolutions of neutrino-heating time scale (left) and volume-integrated net neutrino- 
heating rate (right) in the 2D (green line) and 3D (red line) models (figures taken from Ref.,^^''' 
reproduced by permissions of the AAS). 
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2.2. 3D fully general relativistic simulations with an approximate neutrino transport 

In addition to the 3D effects mentioned in the previous section, impacts of GR 
on the neutrino-heating mechanism stand out among the biggest open questions in 
the supernova theory. It should be remembered that researchers in the pioneering era 
of supernova simulations tackled GR from very early on/^^^ using a newly derived 
formulation. One year after Colgate & White published their seminal paper^^ , 
Schwartz^^^) reported the first fully GR simulation of stellar collapse to study the 
supernova mechanism, implementing a gray neutrino transport in the ID GR hy- 
drodynamics code*-*. Using the GR Boltzmann equations derived by Lindquist,^^®^ 
Wilson^^^) developed a ID GR-radiation-hydrodynamics code including more real- 
istic (at that time) descriptions of the collisional term than the one adopted by 
Schwartz. "^^^^ ID GR hydrodynamical simulations with the so-called leakage scheme 
for neutrino cooling were also performed to explore hydrodynamical properties up to 
the prompt-shock stagnation.^^^^"^^°) These pioneering studies, albeit with a much 
more simplified neutrino physics than today, did provide a bottom-line of our cur- 
rent understanding of the supernova mechanism (see Bruenn et al.^^^^ for a complete 
list of references for the early GR studies). In the middle of the 1980s, Bruenn^^^^ 
developed a code that coupled ID GR hydrodynamics to the MGFLD transport 
with a relativistic correction of order {v/c) and included the so-called standard set 
of neutrino interactions. By the late 1990s, the ultimate ID simulations, in which 
the GR Boltzmann transport is coupled to ID GR hydrodynamics, became feasi- 

l3lg_12),14),121),123)-126) 

Bruenn et al.^^^^ demonstrated clearly that the average neutrino energies of all 
neutrino flavors are higher in GR than in Newtonian gravity during the shock-heating 
phase. They also pointed out that the redshift and gravitational time dilation, the 
agents to counteract, are rather minor. Employing the best weak interaction rates 
available at present, Lentz et al.^^^^ reported very recently the update of Bruenn et 
al.,^^^^ in which they showed that the neglect of the observer corrections in the trans- 
port equation particularly does harm to neutrino-driven explosions. In these full- 
fledged ID simulations, a commonly observed disadvantageous aspect of GR is that 
the residency time of material in the gain region is shorter owing to stronger grav- 
itational pull. All these effects taken into account, GR is negative in the neutrino- 
heating mechanism in ID. In fact, switching from Newtonian to GR hydrodynamics, 
we find that the maximum shock radius becomes ~ 20% smaller in the postbounce 
phase (e.g., 127)). 

In the most advanced multi-D simulations with spectral neutrino transport men- 
tioned earlier, GR effects are addressed at best by a modified gravitational potential 
that is adopted from the ID post-Newtonian correction.^^^'^^^'^^^'^^^^ A possible 
drawback of this prescription is that the total-energy conservation is compromised 



*' Cited from his paper, "In this calculation, the neutrino luminosity of the core is found to be 
10^* erg/s, or 1/2 a solar rest mass per second !! .... This is the mechanism which the supernova 
explodes". The neutrino luminosity rarely becomes so high in modern simulations, but it is surprising 
that the potential impact of GR on the neutrino-heating mechanism was already indicated in the 
very first GR simulation. 



14 



K. Kotake, K. Sumiyoshi, S. Yamada et al. 



owing to the term added artificially to the Poisson equation for self-gravity. Since 
the supernova engine is powered by the gravitational energy, we have to avoid any 
potential inaccuracies in energy conservation. On the other hand, there have been 
a number of fully general relativistic simulations of massive star collapse thus far 
both in 2D (e.g., 129)) and in 3D (e.g., 130), 131), and references therein). The 
so-called conformal-fiatness approximation (CFC) has been also employed. ^^^^'^^^^ 
In these computations the treatment of neutrino transport was overly compromised, 
e.g., with the transfer being entirely replaced by a prescribed Yg formula^^^^ or the 
so-called leakage scheme being adopted.^^^^'-*-^^)*) 

In this section, we present results from our first generation muIti-D hydrodynam- 
ical simulations in full GR that incorporate an approximate neutrino transport. ^^^^ 
The code is a marriage of an adaptive- mesh-refinemcnt (AMR), conservative 3D 
GR MHD code developed by Kuroda and Umeda,^^'^-' and the approximate neutrino 
transport code that we newly developed in this work. Our GR code is based on the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism. ^'^^-''^^^^ Hydrodynamics 
can be solved cither in full GR or in special relativity (SR) , a feature that allows us 
to investigate purely GR effects on the supernova dynamics. Using the so-called Ml 
closure scheme with an analytic variable Eddington factor, we solve the radiation 
energy and momentum. This part of the code is partially based on the Thome's 
momentum formalism, which was recently extended by Shibata et al.^^^^ so that 
it should be more suitable for neutrino transport. To simplify the source terms of 
the transport equations, on the other hand, a multi-flavour neutrino leakage scheme 
is also employed partially. The new code is designed so that it could evolve the Ein- 
stein equation and GR radiation-hydrodynamical equations self-consistently, satisfy- 
ing the Hamiltonian and momentum constraints. The AMR technique implemented 
in the 3D code enables us to follow the dynamics from the onset of gravitational 
core-collapse of a 15 Mq star through bounce up to ^ 100ms after bounce in this 
study. We compute four models with different combinations of SR/GR and 1D/3D, 
which we label as ID-SR, ID-GR, 3D-SR and 3D-GR, respectively. Limited to the 
early postbounce phase (t < 100ms), we discuss exploratory results in the following 
sections to illuminate GR effects in the multi-D neutrino-heating mechanism. 

2.2.1. Hydrodynamical features in full GR 3D simulations 

Four snapshots in Fig. 7 are helpful to characterize the postbounce features in 
our 3D-GR model**^ . The top left panel shows the distribution of entropy per baryon 
at t ~ 10ms, when the bounce shock stalls at a radius of ~ 90km (shown as a central 
blueish sphere). Comparing the top left with top right panel in Fig. 7, we see that the 
shock (a greenish sphere in the top right panel) becomes bigger. This implies that 
the bounce shock turns into a so-called "passive" shock, which expands outward 



*-* Very recently, Miiller et al.^'"*^' reported explosions for 11.2 and 15Mq stars based on their 
2D GR simulations in CFC with detailed neutrino transport similar to 36) being implemented.^^*-' 
**' The 3D computational domain is a cube of 10000^ km^ volume fit in the Cartesian coordinates. 
The maximum refinement level in AMR is 5 at the beginning and then incremented as the collapse 
proceeds. The criterion for incrementation is renewed when the central density exceeds IQ^'^^^^'^^-^ 
g/cm^ , yielding an effective resolution of Ax ~ 600m at bounce (see 140) for more details) . 
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Fig. 7. 3D snapshots of entropy per baryon at four different times (top left; t — 10ms, top right; 
t — 40ms, bottom left; t — 80ms, and bottom right; t = 100ms) for model 3D-GR. The contours 
in the a; = 0, y = 0, and z — planes are projected on the back right, back bottom and back 
left sidewalls, respectively, to visualize the 3D structures. In each plot, an arbitrarily chosen 
iso-entropy surface is displayed. The linear scale is indicated along each axis in unit of km 
(figures taken from Ref.139), reproduced by permissions of the AAS). 

gradually with all matter advecting inward after passing through the shock (e.g., 
36)). At this stage, there forms a gain region, in which neutrino- heating dominates 
over local cooling. The neutrino-driven convection gradually develops from this 
point on. The sidewalls in the top right panel also demonstrate the growth of the 
postshock convection. The entropy behind the standing shock becomes higher with 
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Fig. 8. Evolutions of average shock radii as a function of time after bounce for models ID-SR (green 
line), 3D-SR(blue line), ID-GR (black line), and 3D-GR (red line), respectively. The maximum 
shock radii and the time scale of shock recession for model ID-GR is similar to those obtained in 
the previous ID simulations for the same progenitor model that incorporated neutrino transport 
more in detail (e.g., 127), 138)). 

time owing to the neutrino-heating, which can be inferred from yellowish bubbles in 
the bottom left panel. These high entropy bubbles {s > lOke) rise and sink behind 
the standing shock. The shock deformation is dominated by unipolar and bipolar 
modes, which may be a characteristic feature of the SASI. The neutrino-heated region 
becomes larger with time in a non-axisymmetric way, which is evident in the bubbly 
structures that are shown as reddish regions in the bottom right panel. 

In these simulations up to 100 ms after bounce, the largest shock radius is 
recorded in model 3D-GR (red line in Fig. 8). The other models have already seen the 
shock recession by this time. Before we focus on the reason for it in the next section, 
let us compare the activities of convection and SASI between the computations in SR 
and GR. Figure 9 displays the angle-averaged Brunt- Vaisala (BV) frequency u-Qy^^^ 
(left panels) and pressure dispersion Ap*^ in a logarithmic scale (right panels) for 
models 3D-SR (top two panels) and 3D-GR (middle two panels), respectively. 

Irrespective of SR or GR, there are typically three convectively unstable regions 
in the postbounce phase: (1) the greenish region behind the shock**^ at t < 20ms 
that corresponds to the so-called prompt convection, (2) a narrow horizontal strip 
behind the shock that corresponds to the convection sometimes referred to as Bethe 
convection, (3) a thick horizontal strip above the PNS at a radius of r ~ 10 — 

*' This is defined as Ap = — — , where (^4) represents the angle average of quantity A. 

Note that the shock is indicated by a white thin line that rises quickly after bounce and 
declines after the passive shock stalls at a radius of r ~ 150km . 
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Fig. 9. Postbounce evolution of angle-averaged Brunt- Vaisala frequencies {ujbv in ms~^) (left pan- 
els) and pressure dispersions Ap (right panels) for models 3D-SR (top row) and 3D-GR (middle 
row), respectively, and the maximum pressure dispersion Zipmax (bottom row). Only convec- 
tively unstable regions (i.e. ujbv > 0) are shown and the white lines represent the boundaries of 
convective regions with lubv = in the left panels. In searching the maximum value of pressure 
dispersions for the bottom panel, we restrict the region to 20 < r < 50km, i.e., the vicinity of 
the coupling radius (figures taken from,^^^' reproduced by permissions of the AAS). 

20km that emerges at t > 60ms. Comparison between the two panels m the left 
column for models 3D-SR and 3D-GR shows that the unshocked core (the central 
part of PNS surrounded by the convective region) is more compact in the GR model 
at t > 50ms. The PNS convection develops only very weakly before t ~ 60ms. 
This is common to both SR and GR cases and due to the stabilizing effect of a 
positive entropy gradient prevailing outside the PNS surface (r ~ 10km). The PNS 
convection becomes vigorous gradually with time afterward as the negative lepton 
gradient develops in the nascent PNS. 
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Next we pay attention to the right two panels of Fig. 9 to infer the activities 
of SASI. In these panels, we may recognize two horizontal strips: one is colored in 
red and shows strong pressure perturbations behind the shock; the other is colored 
in green and roughly corresponds to the bottom of cooling layer that recedes from 
r ~ 80km to r ~ 30km gradually in time from t ~ 30ms to t ~ 100ms. The accreting 
flows that advect from the standing shock on to PNS receive abrupt decelerations 
near the bottom of the cooling layer. Strong pressure perturbations are produced 
there as mentioned above and propagate outward subsequently until they hit the 
shock. The up-going stripes in the figure seem to indicate these outward propagation 
of pressure waves. The features just mentioned may be reminiscent of the so-called 
advectic-acoustic cycle (e.g., 145)-147) and references therein) and are common to 
the SR and GR models. 

The bottom panel of Fig. 9 compares the maximum pressure dispersions that 
the advecting vortices produce in the vicinity of the deceleration region. As is 
evident in the figure, the maximum pressure dispersion is generally larger in the GR 
model (red line) than the SR counterpart (black line) in the early postbounce phase 
(i < 100ms) we study here. This is presumably because stronger gravitation pull in 
GR makes the coupling radius smaller, leading to the production of more energetic 
acoustic waves. Although it is not straightforward to say something very solid only 
from this figure, what we observed so far in our 3D-GR model, i.e., the generation 
of stronger acoustic waves and larger shock radii in GR, suggests at least that 3D is 
not unfavorable for the neutrino-driven explosion. We now move on to more detailed 
discussions on potential impacts of 3D and GR on the neutrino-heating mechanism. 

2.2.2. 3D versus GR: impacts on the neutrino-heating mechanism 

Recalling that the neutrino heating rate can be expressed as cx Lj^(e^),^^^) 
we first analyze the neutrino luminosities (L^,) and mean energies {{Eh)). We then 
compare the dwell time with the neutrino-heating time in the gain region and discuss 
which one, 3D-SR or 3D-GR, is more likely to satisfy the criterion for shock revival. 




Fig. 10. Luminosities of all neutrino flavors for all the models as a function of time. Ve and 
are displayed in the left panel and is presented in the right panel (figures taken from,^^®' 
reproduced by permissions of the A AS). 



*' Scheck et al. ' referred to this region as the "coupling radius", at which the coupling of 
vortices and acoustic waves takes place. 
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Figure 10 shows for all the computed models the evolutions of neutrino lumi- 
nosities of all the species: Ve, Vx (left panel) and Ve (right panel). The spike in 
the Ve luminosity is the so-called ncutronization burst that occurs when the shock 
passes through the neutrino sphere for i/g. The peak luminosities for the GR 
models are L^,^ ~ 3 x lO^^erg s~^, slightly larger than those for the SR models, 
{Lij^ ~ 2.9 X lO^^erg s"-*^) but are rather insensitive to dimensionality. This trend is 
qualitatively similar to what was found in 121). On the other hand, recent studies 
with weak interactions being treated in more detail in approximate Boltzmann trans- 
port have demonstrated that the peak luminosity is ~ 10% smaller in GR than in 
Newtonian gravity (e.g., 127), 138)). This may carry an important message that the 
Boltzmann transport should be implemented in the full GR simulations to obtain a 
~ 10% accuracy, which is not small at all when speaking about the neutrino-heating 
mechanism. 

After the neutronization burst {t > 10ms), the luminosities increase with 
time in the GR models whereas they arc almost constant in the SR models in the 
early postbounce phase up to t ~ 100ms (green and blue lines in Fig. 10). The 
luminosities are highest in model 3D-GR after t ~ 50ms (red line in the right panel 
of Fig. 10). This is also the case for the luminosities (left panel). Although the 
luminosities change with time, the luminosities of different neutrino flavors satisfy 
the following orders in general: 

ue-. 3D-GR > IDGR, 3D-SR ~ ID-SR, 

i^e- 3D-GR > IDGR, 3D-SR > ID-SR, 

u^: 3D-GR > IDGR, 3D-SR > ID-SR. 
In short, both 3D and GR raise the neutrino luminosity in the early postbounce 
phase. More specifically, the maximal boost by GR, ~ 50%, is obtained for u^; in 
3D as is found in the left panel of Fig. 10 whereas the maximum gain by 3D is less 
than ^ 20%, which is obtained for Ug in the comparison between models 3D-GR and 
ID-GR. These results indicate that GR holds a comparatively more important key 
to the neutrino luminosity. 

The top two panels in Fig. 11 present the angle averaged RMS neutrino energies 
for fg (left panel) and z^g (right panel) after the ncutronization burst (t > 10ms). 
We obtain the highest energies in model ID-GR (black line) and the second highest 
in model ID-SR. Then comes model 3D-GR followed by model 3D-SR. In accord 
with the previous ID results,^^^-''^^^^'^^^-''^^^^ our 3D models (albeit limited to the 
early postbounce phase) support the expectation that we will obtain higher neutrino 
energies when switching from SR to GR. The deeper gravitational well in GR is the 
reason for the higher neutrino energies. In fact, PNS becomes more compact and, 
as a consequence, hotter in GR, which then leads to smaller and hotter neutrino 
spheres. This is evident when one compares the radii of neutrino spheres between 
the GR and SR models in the bottom panels of Fig. 11. Smaller neutrino energies 
in the 3D models compared with the corresponding ID counter parts (top panels) 
are due to larger neutrino spheres in the former (bottom panels). In fact, the shock 
reaches larger radii in 3D, assisted by the convection and SASI (e.g.. Fig. 8), which 
also helps shift the positions of neutrino spheres outwards. The enlarged neutrino 
spheres in multi-D models are qualitatively consistent with the 2D post-Newtonian 
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results by 36), which incorporated the advanced neutrino transport. 
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Fig. 11. Evolutions of the angle averaged RMS energies (upper panels) and the radii of neutrino 
spheres (lower panels) for (left panels) and Vc. (right panels). The colors of Unes are the same 
as those in Fig. 10 (figures taken from,^''^' reproduced by permissions of the AAS). 

Not ah GR effects are good for the neutrino-heating mechanism. In fact, stronger 
gravitational pull in GR tends to shorten the residency times of accreting matter in 
the gain region. In the following we discuss whether the net GR effect, after all these 
effects being taken into account, is positive or not in the multi-D context. Figure 12 
shows the ratio of the residency time scale to the neutrino-heating time scale for all 
the computed models. Although the computed time ~ 100ms is way too short for 
the stalled shock to be revived, our results clearly suggest that the shock revival is 
most likely to occur in model 3D-GR (red line). Models 3D-SR, ID-SR and ID-GR 
follow it in this order. Thanks to larger degrees of freedom, the residency time scale 
is much longer in the 3D models than in the ID models. In addition, the increases in 
neutrino luminosity and RMS energy via the GR effects (Fig. 11) raise the time scale 
ratio by a factor of < 2 in model 3D-GR (red line) from the SR counter part (blue 
line). Our results hence suggest that the combination of 3D and GR will provide the 
most favorable condition for the neutrino-driven explosion. 

It is expected from Fig. 12 that the shock revival will never occur in the ID 
models, which have already shown the sign of a rapid shock recession by the end 
of the simulations. On the other hand, the time scale ratio remains high in the 
3D models for the last ~ 30ms before the simulations are terminated. For the 
15 Mq progenitor employed in this study, it is expected that the neutrino-driven 
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Fig. 12. The ratio of the residency timescale to the heating timescale for the set of our models as 
functions of post-bounce time (see text for the definition of the timescales). This figure is taken 
from,^^®' reproduced by permissions of the AAS. 

explosions take place at t ~ 200ms at the earliest'^'') and that it might be delayed to 
t > 600ms^^) as already mentioned. The parametric explosion models showed that 
the earlier the shock revival occurs, the stronger the explosion becomes^*^^'^^-'). The 
shock revival times obtained in the previous 2D simulations^^)' ^"^'^^^'^^^ could have 
been shorter if the combination of GR and 3D had been included. We anticipate that 
this can be a possible remedy to turn the relatively under-powered 2D explosions 
into more powerful ones. It is worth pointing out that the combination of GR and 
3D*) should affect not only the supernova dynamics, but also the observational multi- 
messenger signatures (e.g., Ref. 74) for a recent review), such as gravitational- waves 
(e.g., 107), 108), 149)-151)), neutrino emission (e.g., 152)-154)), and nucleosynthetic 
yields (e.g., 155), 156)). To give reliable predictions to these important observables, 
the multi-energy and multi-angle neutrino transport should be incorporated in full 
GR simulations together with more detailed weak interactions. This work is only a 
very first step on the long and winding road. 

§3. Progresses in Boltzmann neutrino transport 

3.1. Numerical simulations with Boltzmann equations: overview 

We briefly overview here the recent progresses in the numerical treatment of 
neutrino transfer with exact Boltzmann equations in the CCSN simulations. ^^^)' ^^^^ 
Although 3D computations of hydrodynamics are now made practicable thanks to 
large computing resources available these days, the neutrino transport in three spa- 
tial dimensions is still a great challenge. It is required to solve the time evolution of 
neutrino distributions in the six dimensional phase space with three components of 

*-* It should be mentioned that MHD effects also remain to be studied (e.g., 63),64),69)-73),140) 
and see also 2) for collective references). 
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neutrino momentum (an energy and two angles) in addition to three spacial dimen- 
sions. Even in spherical (axial) symmetry, three (five) dimensional computations are 
needed. Although various approximations have been proposed so far, solving exact 
Boltzmann equations is highly recommended so that we could remove uncertain- 
ties in the neutrino transfer, the key ingredient in the neutrino-heating mechanism. 
Simplifications such as just dropping energy or angle dependence are not reliable 
in principle, since neutrino interactions are strongly energy-dependent and angular 
distributions are essentially important to accurately estimate neutrino-heating rates 
in the semi-transparent region. 

Under spherical symmetry, the direct solution of the Boltzmann equations for 
neutrino transfer is now possible^^-*' ^^^^ even in GR^^^'^^^' ^^^^ with current 

computing resources. In fact, with these first-principle-based codes, the influences of 
EOS and various neutrino reactions on the supernova dynamics have been examined 
in detail over the years.^^^'^^^'^^^)'^^^^ As mentioned already, it has been consistently 
demonstrated that no explosion is obtained under spherical symmetry. It is well es- 
tablished through the improvements in the numerical treatment of neutrino transport 
during these years, however, that the accurate computation of neutrino transfer is 
indispensable to determine the luminosity and energy spectrum of neutrinos emitted 
from the PNS and the heating rates behind the stalled shock, which will in turn affect 
the shock revival. The Boltzmann neutrino transport is indispensable also for 
a reliable theoretical prediction of neutrinos signals,^^^'^^)'^^)'^^^)'-*^^^) which should 
be compared with future observations by terrestrial neutrino detectors. ^^^'^^^^"-"^^^^ 

With an assumption of axisymmetry, very elaborate, state-of-the-art, approxi- 
mations have been developed in the last couple of years. The flux limited diffusion 
(FLD) method^s)' ^^s). m and the " ray-by-ray" extension of the ID Boltzmann trans- 
port scheme^^)''^^) have been extensively employed to investigate multi-D effects on 
the explosion mechanism. Each approximation has its pros and cons: in the flux- 
limited diffusion approximation, for example, the transport in the semi-transparent 
region is not very reliable; in the "ray-by-ray" approximation, on the other hand, the 
neutrino transfer equations are solved along each radial ray independently and, as 
a consequence, although the computations are highly efficient, the forward-peaked 
distributions of neutrinos in the transparent region tend to be overestimated. Some 
of more recent works combine approximate ID transport schemes such as FLD or 
IDSA with the ray-by-ray technique. ^''-'"^^^ These 2D simulations have demonstrated 
the critical role of hydrodynamical instabilities such as convections and SASI in the 
neutrino-heating mechanism. It should be also mentioned that the exact Boltzmann 
equations were also solved in 2D by 172), 173) with the discrete-ordinate method, 
albeit for a limited number of models. 

Spatially 3D simulations of core-collapse are still in its infancy. In most of 
the earliest 3D simulations^^)' ^^^'^^^'^^^^ the neutrino transfer were just neglected 
or simplified considerably and the authors paid attention to novel features of 3D 
hydrodynamics, in particular instabilities such as convections and SASI. In the so- 
called light bulb approximation, for example, the neutrino luminosity and energy 
spectrum are not solved but prescribed parametrically to seek favorable conditions 
for shock revival. ^^^'^^^^ More recently, by combining the ray-by-ray approximation 
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with FLD^*^) or IDSA,^^-'''^^-''^^^) 3D hydrodynamical simulations with spectral neu- 
trino transport were performed. These 3D computations heralded a new stage of 
supernova simulations. As already mentioned in section 2.1, these models did not 
have sufficient resolution yet and further improvements are needed. In this section 
we will make an attempt to go beyond such approximate neutrino transport schemes 
and solve the exact Boltzmann equations in 3 spatial dimensions, i.e., in 6D phase 
space. For the moment, relativity is neglected in the multi-D Boltzmann transport. 

But before going to the multi-D transport, we will first summarize briefly our ID 
GR radiation-hydrodynamical core-collapse simulations with a Boltzmann solver to 
demonstrate what insights can be obtained into microphysics with these simulations. 
We then report our recent progresses in the coding of multi-D Boltzmann solver, 
presenting the results of some test calculations. 

3.2. ID GR neutrino-radiation hydrodynamics with a Boltzmann solver 

The first-principle simulations by solving the GR hydrodynamical equations and 
exact Boltzmann equations for neutrinos under spherical symmetry enable us to 
avoid uncertainties in numerics and investigate physics, in particular, the influences 
of microphysics on the dynamics in a quantitative manner (see also more recent 
works^^)'"'^^'')'^''^)). This is indeed a good example to see the close connections be- 
tween the latest knowledge on nuclear/particle physics in laboratory and the under- 
standing of astrophysical phenomena. We discuss here how EOS at high densities 
impacts the postbounce dynamics and what information on EOS can be extracted 
from them in return. In the following, we summarize our ID results for ISM© and 
4OM0 stars, paying particular attention to the light curve and energy spectrum of 
neutrinos obtained in the long-term postbounce evolution.^^)'^^)'^^)'^^^)'^^^)"^^^) 

3.2.1. Competing effects of nuclear EOS on the core dynamics 

In the case of the 15Mq star, we computed the evolution up to ~ 1 s after 
bounce^^-* to find out the fate of the stalled shock and see the thermal evolution 
of PNS (see also the long term evolutions by 165)). We adopt two EOS's as in 
the previous sections, i.e., Shen's EOS^^)"^^) and Lattimer&; Swesty's EOS with an 
incomprcssibility of ISOMeV.^^^ We found that neither EOS produced explosions 
and that the shock radii are rather similar in the two cases despite the different 
features of the EOS's (Shen's EOS is harder than Lattimer & Swesty's). It turns 
out that a number of effects are counteracting each other. 

On one hand, the larger symmetry energy of Shen's EOS leads to a smaller 
abundance of free protons, which then reduces electron captures^^^^ and make the 
inner core at bounce more massive, as can be seen in Fig. 13. Note that a larger inner 
core is favorable for explosion. The difference in core mass amounts to ~ O.IM©, 
which can sap the shock energy of ~ 10^^ erg by dissociation of heavy nuclei. On 
the other hand, stiffer Shen's EOS produces core bounce at a slightly lower density 
and, as a consequence, gives a lower core temperature than Lattimer & Swesty's. 
This then reduces the neutrino luminosity and leads to lower heating rates behind 
the shock as demonstrated in Fig. 13. This way the advantage earned during the 
collapsing phase is almost canceled out in the post-bounce phase. Note that the 
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Fig. 13. The velocities at core bounce (left) and heating rates &t t = 150 ms after bounce (right) 
in the collapse of ISM© star. The red and blue lines show the results obtained with Shen's and 
Lattimer & Swesty's EOS's, respectively. 



electron captures on nuclei may dominate over those on free protons. ^^^^ In order 
to address this issue appropriately, we need to implement a multi-nuclei EOS at 
sub-nuclear densities. ^'^•''^^•''^^•''^^''^ Note that the standard EOS's (including Shen 
and Lattimer & Swesty) employ the so-called single nucleus approximation, in which 
an ensemble of heavy nuclei is represented by a single, supposedly most abundant 
nucleus. Note again that the most abundant nucleus need not have the largest 
electron capture rates. In fact, one has to take an ensemble average of the electron 
capture rate multiplied by the abundance of individual nucleus. It is also stressed 
that microphysics such as EOS and weak interaction rates is important also in multi- 
D simulations, since they set the initial shock energy and the emission and absorption 
of neutrinos just in the same way we have seen above. 

3.2.2. Neutrino signals from failed supernovae: extraction of information on EOS 

The long-term simulations of core-collapse of the 40Mq star provides us with an 
opportunity to study black formations and neutrino signals from them. In fact, since 
the size of Fe core is much larger compared with that for ISMq star, it is expected 
that there is no chance of explosion and that continued mass accretion from outer 
envelopes will eventually cause the second collapse to black hole.^^^^^^^^^ Indeed the 
mass of PNS increases rapidly and it reaches the critical mass in ~ 1 s and triggers 
the dynamical collapse to the black hole.^^-*'^^-*'^"^^^ During the thermal evolution 
of the central object, the neutrinos are copiously emitted by electron and positron 
captures as well as thermal productions. The neutrino emission is terminated soon 
after the event horizon is formed and the emission region is swallowed into it. Hence 
the duration of the neutrino emission is essentially determined by the time of the 
second collapse of the accreting PNS into the black hole. 

In order to specify common and different characters of hydrodynamics and neu- 
trino emission in the black-hole forming collapse, we investigated other progenitor 
models in the mass range of 40— SOMq.^^^'^^^^'^'^^-' In addition to the standard EOS's, 
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we also employed hyperonic EOS^'''^^' ^"^^^ as well as quark EOS.^^^ In Fig. 14 we show 
the comparison of the three EOS's: Shen's EOS, Lattimer & Swesty's EOS (with 
an incompressibility of ISOMeV) and hyperonic EOS. It is evident that the energies 
and luminosities of neutrinos rise rapidly due to the increase of temperatures inside 
the slowly contracting PNS and the persistent mass accretion. The duration of the 
neutrino emission is only 0.6-1.3 s. These features are different from those for the 
ordinary neutrino emission in supernovae, which last ~ 20 s with gradually decreas- 
ing energies and luminosities. Hence, it is possible to distinguish the black hole 
formations from the neutron star formations by the neutrino signals. In fact, taking 
properly into account the detector properties as well as neutrino oscillations, we es- 
timate that Super-Kamiokande will record ~10^ events for a galactic event, 
which are comparable to those for ordinary supernovae. ^^^^ 




Fig. 14. Time profiles of the neutrino emission in the collapse of the 4OM0 star. The neutrino 
luminosities (left) and average energies (right) of three neutrino species are shown for three 
EOS's: Shen's EOS, Lattimer & Swesty's EOS and Hyperonic EOS. 



As is evident in Fig. 14, the features of neutrino emission in the black-hole form- 
ing collapse are sensitive to EOS at supernuclear densities. A softer EOS ends up 
with shorter neutrino emission, since the critical mass for the second collapse to black 
hole is lower. ^^^'^^^^ The reason for the softening of EOS may be due to interactions 
between baryons or to the emergence of new degrees of freedom such as hyper- 
Q]^gi77),i78) Qj, quarks. We further analyzed the neutrino signals obtained with the 
rather soft nucleonic EOS (Lattimer &: Swesty'S EOS with an incompressibility of 
220MeV) and with the Hyperonic EOS^^^ by a statistical method. We demonstrated 
that although the durations of neutrino emission are similar to each other, they are 
still distinguishable by Super-Kamiokande if they occur in the Galaxy. These differ- 
ences of signals in turn can be used as a useful probe into the EOS at very high den- 
sities, which may be accessible to the next-generation terrestrial experiments. ^^^-^ 
It should be mentioned that the features of dynamics and neutrino emission depend 
also on the profile of progenitors mainly through the accretion rate.^''^'^''^^ Hence, 
it is important to perform more systematic simulations and construct the templates 
of signals that can be compared with future observations.^^''''' "^^'^^'^^^^'^^^^ 
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3.2.3. Parallel computing in ID simulations 

It is profitable to comment on the aspect of supercomputing in the ID studies. 
In our numerical simulations described above, the numerical inversion of block tridi- 
agonal matrices that appear in the discretization and linearization of basic equations 
is the major computational load. We implemented the block cyclic reduction^^^^ as 
an efficient parallelled algorithm of matrix inversion that replaces the conventional, 
serial algorithm of the Feautrier method. ^^^^ The numerical technique developed for 
ID spherically symmetric neutrino transfer can be easily extended to multi-D by the 
use of the approximate ray-by-ray implementation of Boltzmann solvers. 

3.3. Multi-D neutrino transport: numerical solutions of 6D Boltzmann equations 

Thanks to recent expansions of supercomputing resources, it has now become 
feasible to numerically solve the Boltzmann equations for neutrino transfer in three 
spacial dimensions. We have indeed developed a Boltzmann solver with multi-energy 
and multi-angle groups that is meant for multi-D simulations. ^^^^ The code is based 
on the so-called discrete-ordinate (S„) method in six dimensions (see Ref. 189) for 
an alternative approach by Monte Carlo scheme). 

We describe the neutrino distribution in the space coordinate with radial N^- 
, polar Ng-, and azimuthal N(f,-gnd points and in the neutrino momentum space 
with energy A^^-grid points and angle Ng^- and A^^^-grid points. A fully implicit 
differencing is adopted for time advancement. We choose the inertial frame to write 
down the Boltzmann equations, in which the advection terms have the simplest 
expressions. Then we need to consider the Lorentz transformations to evaluate the 
collision terms, which become simplest in the comoving frame. For the moment, 
however, neglecting all the corrections of the order of v/c or higher, we do not 
distinguish these two frames*^. 

The basic set of neutrino reactions"'^'^-''"^^^^ including pair processes but not in- 
elastic scatterings is implemented in the collision terms. Three species of neutrinos 
(fe, i>e, ^^/t) treated. The standard EOS tables are employed to obtain the ther- 
modynamical quantities and composition of matter, which arc necessary to evaluate 
the collision terms. In the following test computations we adopt Shen's or Lattimer 
& Swesty's EOS. 

3.3.1. Some basic tests 

A suit of numerical tests have been done to validate the newly developed Boltz- 
mann solver. We have first computed multi-D transport in uniform matter both in 
the diffusion and free streaming regimes. In the left panel of Fig. 15, the results of 
2D/3D diffusion of Gaussian packets are displayed. We have found a satisfactory 
agreement with the exact solutions. In the right panel of the same figure, on the 
hand, we have shown the free propagation of neutrinos in a designated direction. 
Substantial numerical diffusion are apparent in this case. Note that this test is too 
demanding and no such situation becomes important in supernova simulations. The 
intermediate regime, the most important one, have been examined by utilizing the 

*^ We have an idea to rigorously treat the Doppler effects and angular aberrations in the coUi- 
sional integrations and will publish it elsewhere 
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Fig. 15. Examples from the test suit of the time evolution of neutrino distributions by the 3D 
Boltzmann solver. Snapshot of neutrino density from the time evolution of the diffusion of 2D 
Gaussian packet by surface plot (left) and the neutrino beam injected from the point source 
in 3D box (right). The color expresses the neutrino density. The number of grid points is 
NrXNg = 100 X 96 (left) and iV^ x A^e x Af^ = 50 x 48 x 48 (left) with iVs„ x TV^^ x Af^ = 12 x 12 x 4. 



formal solutions for more realistic matter distributions as in the following. 

We prepare artificially deformed spheroidal/ellipsoidal cores based on the post- 
bounce profile obtained in the ID GR neutrino-radiation hydrodynamical core- 
collapse simulations of the ISM© star.^^^ With the background profiles of density, 
proton fraction and temperature being fixed, we follow the time evolution from a 
certain initial condition for a sufficiently long period (-^10 ms) to obtain the steady 
neutrino distributions. In this computation we adopt the Shen's EOS table for the 
evaluation of the collision terms. The neutrino densities and fiuxes are evaluated by 
appropriate integrations of the neutrino distribution functions over the momentum 
space. Various angle- and energy moments of the neutrino distribution functions in- 
cluding the fiux factor and Eddington tensors are also examined (see 158) for detailed 
analyses). Moreover, we have also obtained from the collision terms the detailed in- 
formation on neutrino reactions such as mean free paths, deleptonization rates and 
cooling/heating rates. 

In Fig. 16, we show the results for the axially symmetric, spheroidal core. An 
oblate PNS sits at the center and the shock is standing around 140 — 200km in 
this model. The electron-type neutrinos (top left panel) and anti-neutrinos (top 
right panel) are abundant at and off center, respectively, reflecting the degeneracy 
of electrons. It is as expected intuitively for the oblate shape of the core that the 
radial flux is enhanced near the polar axis. The polar fluxes (^-component of the 
flux vectors) are substantial at intermediate polar angles, since the neutrino fluxes 
are inclined towards the polar axis. 

In Fig. 17, we show the 3D case, in which the core is deformed to a non- 
axisymmetric, ellipsoidal shape by adding an azimuthal {(f)) dependence to the spheroidal 
conflguration employed above. Displayed in the flgure are densities and fluxes of 
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Fig. 16. Color maps of neutrino densities and fluxes in tfie lialf of meridian slice for an axisymmet- 
rically deformed spheroidal core. The top panels display the densities of electron-type neutrinos 
(left panel) and anti-neutrinos (right panel) whereas the bottom left and right panels present 
the radial and polar components of flux vectors of electron-type anti-neutrinos, respectively. 
The number of grid points is Nr x Ne x = 200 x 18 x 9 with Ne^ x iV^,^ x A^^ = 6 x 12 x 14. 



three species of neutrinos in the first octant of a meridian sHce with an azimuthal 
angle of (/> = 0.44 radian. The asymmetry between the pole and the equator in this 
slice is slightly smaller than that in Fig. 16. The density of electron-type neutri- 
nos is high at center whereas electron-type anti-neutrinos and mu-type neutrinos 
are abundant off center, where the temperatures are higher than at center. This 
feature is similar to what we saw in the previous axisymmetric case. The fluxes of 
electron-type anti-neutrinos reflect the 3D deformation, though (lower panels of the 
figure). The radial flux around the polar axis is larger than that near the equatorial 
plane although the asymmetry is less remarkable compared with the 2D case. The 
polar fluxes are again non-negligible and largest at intermediate polar angles. Owing 
to the azimuthal dependence of deformation, the azimuthal fluxes are non-vanishing 
and substantial indeed in the broad region around r = 50 — 100km. These non-radial 
fluxes are important to accurately describe the global behavior of neutrino transfer 
in 3D and may affect the neutrino heating rates and, as a consequence, explosions. 
We remark that these non-radial fluxes can be automatically and properly treated by 
the 3D Boltzmann solver unlike the ray-by-ray approximation, in which the neutrino 
fluxes are assumed to be radial. It is also pointed out that the non-radial fluxes are 
non-negligible up to ~ 100km and it is dubious that FLD can give the flux vectors 
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Fig. 17. Color maps of neutrino densities and fluxes in the quadrant of meridian slice for a non- 
axisymmetrically deformed ellipsoidal core. The densities of electron-type neutrinos (left panel), 
anti-neutrinos (middle panel) and mu-type neutrinos (right panel) are shown on the top rows 
whereas the radial (left panel), polar (middle panel) and azimuthal (right panel) components 
of the flux vectors of electron-type anti-neutrinos are shown, respectively, on the bottom row. 
The number of grid points is iVr x iVs x A^^ = 200 x 9 x 9 with Ne^ x Af^„ x iV^ = 6 x 12 x 14. 



correctly there (see also Ref. 172)). 

3.3.2. Several demonstrations for more realistic backgrounds 

We proceed to some more test computations done under more realistic settings. 
We employ the snapshot at 100 ms after bounce that is obtained in the simulations 
of the II.2M0 star^^'^) discussed in §2. In the left panel of Fig. 18, we show the 3D 
entropy distribution in the supernova core. It is clear that the matter distribution is 
deformed both globally and locally due to the hydrodynamical instabilities below the 
shock (shown as a greenish sphere), which is located around 200—300 km. Adopting 
Lattimer & Swesty's EOS for this test and fixing the background, we follow the time 
evolutions of neutrino distribution functions in the 6D phase space from an almost 
vanishing population until time-independent solutions are obtained. The Boltzmann 
solver treats the building up of equilibrium distributions in the optically thick region, 
neutrino cooling and heating in the intermediate region and outward free streaming 
in the optically thin outer layers simultaneously. A snapshot of some surfaces with 
a constant density of anti-neutrinos is shown in the right panel of Fig. 18. It is 
seen that the neutrino distribution is rather spherical near the center and becomes 
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asymmetric globally and locally in the outer layers, tracing the matter distributions 
(compare with the left panel). 




Fig. 18. Iso-entropy surfaces at 100 ms after bounce in the collapse of the II.2M0 star (left panel) 
and surfaces of a constant density of electron-type anti-neutrinos in the evolution by the Boltz- 
mann equations (right). The number of grid points is A'^^ x Ne x = 128 x 16 x 32 with 
Ne^ x 7V^„ x iV. = 6 X 12 x 14. 

Fig. 19 shows the computed density of electron-type anti-neutrinos (left panel) 
and net heating rates (right panel) on a slightly off center slice, respectively. The 
electron-type anti-neutrinos are abundant again in the off center region, where the 
temperature is high and all flavors of neutrinos are produced thermally. The iso- 
density surfaces are prolate, reflecting the matter distributions. The gain radius is 
located around r ~ 100 km and the global asymmetry of heating/cooling regions is 
moderate. The heating/cooling rates are calculated directly from the integrations of 
the collision terms. For implementing the radiation module in a hydrodynamic code, 
these quantities describing the changes in energy and compositions of matter enter in 
the right-hand-side of the hydrodynamic equations (with negative sign) . This is what 
we are currently undertaking (Nagakura et al. in preparation ^^^^). And then, the 
next important task is to implement the velocity dependent terms in the transport 
equations, which is indispensable for the accurate treatment of the relativistic effects. 
After that, we plan to study hydrodynamical instabilities in the supernova core, then 
move on to perform a full-scale simulation starting from the onset of gravitational 
collapse, through core bounce to shock-stall, until shock revival and explosion in 
a consistent manner. In doing so, highly competitive supercomputing resources in 
Japan, in particular the "K computer" *\ the fastest one in the world as of November 
2011, will be helpful. 



*' it is named after the Japanese word of "Kei" , meaning 10 quadrillion (10 petaflops). Note that 
the given name of the corresponding author of this paper has nothing to do with the supercomputer 
(except for the fact that the person is going to use it for SN simulations). 
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Fig. 19. The density of electron-type anti-neutrinos (left) and net heating rate (right) at 100 ms 
after bounce on a slightly off-center slice. The background matter distribution is obtained from 
the 3D simulation of 11.2Mq progenitor and fixed in the computations of neutrino transport. 



3.3.3. Neutrino transport in highly non-spherical environments: collapsars 

As an additional demonstration of the capability of our new 3D code, we take 
the collapsar model for gamma ray bursts, which are more energetic and 
asymmetric than CCSNe. Gravitational collapse of very rapidly rotating, massive 
cores results in the formation of black hole with a surrounding disk, which is presum- 
ably responsible for jet formations. To pin down the physical processes to form the 
relativistic jets in the collapsar model is a long-standing issue. Annihilations of neu- 
trino pairs emitted from the disks are one of the plausible mechanisms (e.g.,^^^^"^^^-* 
and references therein). Mass ejections from the disk and/or jet through neutrino 
interactions are attracting broad attention in the studies of nucleosynthesis of heavy 
elements. ^^^^"^"^"^ Neutrino transport in the collapsar model is hence an indispens- 
able ingredient for the investigations both of jet formations and nucleosynthesis. 
Since the black hole and disk system is highly non-spherical, numerical approaches 
are almost mandatory. Although extensive studies mainly by utilizing a ray-tracing 
technique have been reported so far to this end,^''^)"^''^) the spectral treatment of 
neutrino transfer has been a major undertaking. As a very first step toward bet- 
ter description of neutrino transport in collapsar, we employ our 3D code here to 
obtain time-independent neutrino distributions (see also Ref. 206)) for a matter 
configuration extracted from a hydrodynamical simulation of collapsar. 

Distributions of density, temperature and electron fraction are provided by 
2D GR simulations of gravitational collapse of a IOOMq star with a rapid rota- 
tion. The density profile after black hole formation is shown in the left panel 
of Fig. 20. The black hole sits at the center, which is removed in the neutrino cal- 
culation and shown as a white circle in the figure, and a dense disk surrounds it. 
Not shown explicitly in the figure, there are accretions and outfiows outside the disk. 
Fixing the matter distribution, we compute the time evolutions of neutrino distri- 
bution functions until the steady state is reached. We adopt the Shen's EOS table 



32 K. Kotake, K. Sumiyoshi, S. Yamada et al. 




" Icf"! I R [cm] 



Fig. 20. The density profile in the quadrant of a meridian section used for the computation of 
neutrino transport with the Boltzmann solver (left) and the resultant steady distributions of 
electron-type anti-neutrinos (right). The neutrino densities are shown by a color map in log 
scale and fluxes are presented by arrows in the right panel. The number of grid points is 
Nr X Ne X = 100 x 45 x 3 with Ng^ x N^^ x A^^ = 6 x 12 x 14. 

for this simulation. 

The density of electron-type anti-neutrinos is shown in the right panel of Fig. 20. 
It is found that the electron-type anti-neutrinos are abundant in the outer part of the 
disk, where the temperatures are high, whereas the density of electron-type neutrinos 
(not shown in the figure) is high near the the black hole. The neutrino fluxes (shown 
as arrows in the figure) reflect the geometry of the system, very roughly agreeing with 
the local gradient of neutrino density, and completely non-radial. This is a situation 
that the ray-by-ray approach is certainly inappropriate. Detailed information such 
as the location of neutrino spheres and energy spectra of all species of neutrinos 
is essential to investigate the dynamics of outflows and nucleosynthesis inside them. 
We will address these issues in the near future with the new Boltzmann solver, which 
will be combined with a 3D GR hydrodynamics code described in section 2.2. 

3.3.4. Progresses in the new algorithms for large-matrix inversion 

It is worth mentioning that the exact 3D Boltzmann solver demonstrated above is 
a product of our collaboration with computational scientists, who know how to make 
best use of computing resources. They are specialists indeed in the mathematical 
modeling, algorithms and parallel computing. In this section, we briefly describe our 
recent progresses in this aspect of the development of the 3D Boltzmann solver. 

As mentioned earlier, our scheme is based on the finite differencing of the Boltz- 
mann equations {Sn method), which is implicit in time. The resulting equations are 
a linear system with a large sparse matrix. The main computational load comes from 
the inversion of this matrix, which has to be done at every time step. In the left 
panel of Fig. 21, we show the positions of non-zero elements in the matrix we obtain 
from the discretization of the original equations. It is found that the matrix consists 
of dense sub-matrices along the diagonal line (shown as gray boxes in the figure) as 
well as non-zero elements on six off-diagonal lines (labeled as xi, X2 and X3 in the 
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Fig. 21. Left: pattern of the sparse matrix appearing in the hnear system obtained for the imphcit 
discretization of the Boltzmann equations. N and M denote the numbers of spatial grids (N,. Ng 
N^) and neutrino grids (Na^N0^Ne), respectively. For the studies on the current supercomputers 
without energy couplings, the size of diagonal black matrices (gray) is N^^N^^. Right: number 
of iterations as a function of the time step for different pre-conditioners, i.e., the point Jacobi 
method (blue crosses) and newly developed method (red crosses). The number of grid points 
for the numerical experiment is Nr x Ng x N,i, = 200 x 9 x 9 with Ne^ x N^^ x = 6 x 12 x 14. 

figure). The sub-matrices originate from the collision terms that express local emis- 
sion and absorption of neutrinos as well as coupling by neutrino scattering moving in 
different directions, whereas the off-diagonal lines of non-zero elements corresponds 
to the spatial advection terms. The total number of grid points in the 6D phase 
space amounts to ~ 10^ in a typical simulation. The size of dense sub-matrix is 
~ 100 in the current studies and will be larger for higher angular resolutions and/or 
full energy-couplings. 

For the inversion of matrices of this size, iterative methods'^''^^ are the first choice. 
We employ as a standard option the Bi-CGSTAB method, utilizing a program in the 
Templates library, together with the point-Jacobi method as a pre-conditioner. 
We obtain convergence typically within 20 iterations with a residual error of 10~^. 
This is of course a function of the time step. At. As we increase it, convergence 
becomes slower because the diagonal elements become less dominant. As a matter 
of fact, sometimes no convergence is obtained even after 200 iterations in the simu- 
lations for realistic matter distributions extracted from dynamical models in §3.3.2. 
We certainly need to find a better way to improve convergence. 

Recently we have found a new method to optimize the pre-conditioning. Tak- 
ing out a set of the matrices from our simulations of 3D neutrino transfer, which 
presents the slow-convergence problem in the standard approach with the point- 
Jacobi method, their properties have been analyzed in detail to propose a parameter- 
optimized damped Jacobi-type pre-conditioner, details of which will be published 
elsewhere. ^^'^^ The convergence efficiency is compared between the two pre-conditioners 
for the same matrices extracted from the 3D Boltzmann simulations. In the right 
panel of Fig. 21, we show the numbers of iteration as a function of the time step. At, 
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for a representative case. As mentioned already, the convergence becomes very slow 
for At > 10~^s and no convergence is obtained for At > 3 x 10~^s even after 200 
iterations when the point-Jacobi method is employed. On the contrary, with the new 
pre-conditioning method, the convergence is improved drastically. The time steps 
can be increased by a factor of 100 up to > 10~^s. This is favorable particularly for 
long-term computations. It is true that the computational cost of the new method 
is higher but it is just by a factor of ~ 10 compared with the standard method. Our 
efforts have hence paid off and we have achieved a speed-up by a factor of ^ 10. We 
are currently applying the new method to various cases to see if such a good perfor- 
mance is retained or not. We will also continue to seek for even better methods, since 
we expect that the matrix will be larger in the productive runs of neutrino-radiation 
hydro dynamical simulations in 3D. 

§4. Beyond the "K computer" 

Rapid growth of the supercomputing capability in Japan these years enables 
us to perform large-scale simulations such as those presented above. 3D supernova 
simulations with sufficient resolutions definitely require the K computer and more 

even beyond Exa-flops scale platforms. We remark also that allocations of sufficiently 
long cpu-time on such facilities are also indispensable for long-term computations 
such as those of delayed neutrino-driven explosions. 

As reported in §2 and §3.3.2, the first 3D simulations of corc-collapsc supcrnovac 
with spectral neutrino transport by the ray-by-ray IDSA were performed on the 
currently available supercomputers. It was demonstrated that the numerical grid 
deployed in the computations was not fine enough to draw a solid conclusion on the 
3D explosion mechanism. Scaled-up simulations are scheduled on the K computer 
in Kobe, Japan. The 3D neutrino transfer with the Boltzmann solver requires even 
larger memory and speed specifications. Summarized in Table I are the target sizes 
of numerical grids we deploy for the 3D Boltzmann transport as well as the required 
memory sizes and expected operation numbers per time-step on both the current 
and future platforms. Note that we have to follow more than 10^ time-steps for a 
productive run. 

We assume in the table that the computational load mainly comes from the 
inversion of the sub-matrices that account for the local emission and absorption as 
well as scattering of neutrinos; the operation numbers for the currently adopted 
inversion scheme, i.e., the Bi-CGSTAB method with the preconditioner discussed 
above, are proportional to N^N^^^i^, in which and N angle are the number of energy 
and angle grid points for neutrino transport, respectively; the latter gives the size 
of the dense sub-matrices on the diagonal line, since the scatterings couple different 
angular grid points. It is found that on the currently available supercomputers, we 
can afford only moderate resolutions for the 3D transfer, limited by the necessary 
memory size, which amounts to 2TB for the storage of the sparse matrices in the 
linear system (plus 20GB for the distribution functions of three species of neutrinos) . 
The number of floating-point operations per time-step is estimated to be 6 x 10^^ 
for a single species of neutrinos. 
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Table I. Target grid sizes and required computational resources for the current and future su- 
percomputers. The currently available (typical) supercomputers, the K computer, Exa-scale 
supercomputers and beyond are listed. The numbers of radial, polar and azimuthal grid points 
in the space (Nr, Ne, N^) and those of angles and energy grid points in the momentum space 
(Ne„, N^^, Ne) are given together with the required memory sizes for the storage of the dis- 
tribution functions of 3 species of neutrinos (4th column denoted as fv's) and the matrix in 
the system of linear equations (5th column denoted as Matrix) as well as the floating-point 
operations (6th column denoted as Operations). The last row refers to the case, in which energy 
couplings by inelastic scatterings are fully taken into account. 



Platforms 


Space (N^NoN^) 


Neutrino (Ne^N,^„Ne) 




Matrix 


Operations 


Current 


256 X 32 X 64 


8 X 12 X 14 


2 X 10' GB 


2 TB 


6 X 10'^ 


K computer 


512 X 64 X 128 


12 X 24 X 20 


6 X 10^ GB 


2 X 10^ TB 


2 X 10'^ 


Exarscale 


512 X 128 X 256 


24 X 24 X 24 


6 TB 


3 PB 


8 X 10^® 


Ei,-coupling 


512 X 128 X 256 


24 X 24 X 24 


6 TB 


8 X 10^ PB 


4 X 10^^ 



With the K computer, which is now in operation, we can deploy twice finer 
spatial, energy and angle grids; the memory size is 2 x 10^ TB for the storage of 
matrices and 6 x 10^ GB for the neutrino distribution functions; the floating-point 
operations per time-step becomes 2x10^^, which certainly requires lOPflops-class 
supercomputers. Even with the K computer, though, it will be hard to follow the 
evolution of supernova cores over 1 s in 3D and such long-term simulations may be 
limited to axisymmetric 2D models. Nevertheless, the K computer makes it possible 
for us to accomplish systematic, high-resolution 2D Boltzmann simulations, which 
are still important on their own right. 

To perform long-term 3D simulations by the 3D Boltzmann solver, on the other 
hand, wc need supercomputers of Exaflops scale, which are expected to come as a 
next-generation platform. As shown in the table, the 3D neutrino radiation hydro- 
dynamics simulations with a sufficient resolution will be feasible only if a memory 
of 3 PB is available and the computational speed is fast enough to handle 8 x 10^^ 
operations per time-step. So far we have ignored inelastic scatterings, which would 
couple different energy grid points and increase the size of the dense sub-matrices 
by a factor of N^. Different energies are also coupled by Doppler effect and gravita- 
tional redshift, which are neglected in the current version of our Boltzmann solver. 
If these effects are taken into account and the resultant enlarged sub-matrices are 
to be inverted in the same way, the required memory and operation numbers are 
gigantic as given in the last row of the table, since they are proportional to and 
iV|, respectively. Last but not least, the implementation of GR in the Boltzmann 
solver, the ultimate goal of our project, should be addressed at an appropriate point 
during this scale-up. 

It is now obvious to readers that the supernova research is a subject of su- 
percomputing science that keeps step with the advancement of hard and softwares 
for supercomputing. Hopefully, the next generation supercomputers will provide us 
with the opportunity to finally reach the goal. We hope also that our quest for the 
supernova mechanism will in turn contribute to pushing the limit of computational 
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science in the decade to come. 
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